{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Boson Sampling with MPS\n",
"\n",
"In this notebook, we explain how to use the MPS (Matrix Product State) backend to simulate a linear circuit. MPS simulation is based on a type of tensor network simulation, which gives an approximation of the output states [1][2]. It does the computation on each component of the circuits one-by-one, and not on the whole unitary. The states are represented by tensors, which are then updated at each component. These tensors can be seen as a big set of matrices, and the approximation is done by chosing the dimension of these matrices, called the *bond dimension*. For this example, we simulate a simple boson sampling problem, with 6 modes and 3 photons [3]."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import perceval as pcvl\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Definition of the circuit\n",
"\n",
"Just like in the Boson Sampling notebook, we generate a Haar-random unitary and its decomposition in a circuit :"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6 modes triangular Boson Sampler :\n"
]
},
{
"data": {
"text/plain": "",
"image/svg+xml": "\n"
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = 6\n",
"unitary = pcvl.Matrix.random_unitary(m)\n",
"mzi = (pcvl.BS() // (0, pcvl.PS(phi=pcvl.Parameter(\"φ_a\")))\n",
" // pcvl.BS() // (1, pcvl.PS(phi=pcvl.Parameter(\"φ_b\"))))\n",
"linear_circuit = pcvl.Circuit.decomposition(unitary, mzi,\n",
" phase_shifter_fn=pcvl.PS,\n",
" shape=\"triangle\")\n",
"\n",
"print(m, \" modes triangular Boson Sampler :\")\n",
"pcvl.pdisplay(linear_circuit, compact = True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## MPS simulation\n",
"\n",
"Let us now define the MPS simulation, using the MPS backend in Perceval."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"simulator_backend = pcvl.BackendFactory().get_backend(\"MPS\")\n",
"simMPS = simulator_backend(linear_circuit)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now choose the size of the matrices (the bond dimension) for our simulation."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"chi = 8\n",
"simMPS.set_cutoff(chi)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And finally, we get the output probability distribution from a given input state with 3 photons."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": "",
"text/html": "
\n\n
state
probability
\n\n\n
|0,0,1,0,2,0>
0.069394
\n
|0,0,2,0,1,0>
0.061578
\n
|1,0,0,1,1,0>
0.059484
\n
|0,1,0,1,0,1>
0.053757
\n
|1,1,1,0,0,0>
0.042388
\n
|0,0,0,2,1,0>
0.039433
\n
|1,0,0,2,0,0>
0.039077
\n
|0,2,0,1,0,0>
0.037225
\n
|0,1,0,2,0,0>
0.031791
\n
|0,1,1,0,1,0>
0.028824
\n
|1,1,0,0,1,0>
0.028661
\n
|0,0,3,0,0,0>
0.028538
\n
|0,1,1,1,0,0>
0.024825
\n
|2,0,0,0,1,0>
0.024031
\n
|0,1,1,0,0,1>
0.022896
\n
|0,1,2,0,0,0>
0.021384
\n
|1,1,0,0,0,1>
0.021065
\n
|0,1,0,0,1,1>
0.020496
\n
|0,2,1,0,0,0>
0.020191
\n
|0,0,0,3,0,0>
0.020153
\n\n
"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n = 3\n",
"input_state = pcvl.BasicState([1]*n + [0]*(m-n))\n",
"\n",
"probs = pcvl.BSDistribution({state : prob for state, prob in simMPS.allstateprob_iterator(input_state)})\n",
"pcvl.pdisplay(probs, max_v=20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Certification of the MPS method :\n",
"\n",
"As we make an approximation by chosing the bond dimension, we have to check when does this approximation becomes good enough. Unfortunately, there is no formula giving the minimal size for a given approximation error. What we can do though is to compute the *Total Variance Distance* (TVD) between an ideal simulation of Boson Sampling, and an approximated one. To compute the ideal one, we can for instance use the *SLOS* backend on Perceval :"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": "",
"text/html": "
\n\n
state
probability
\n\n\n
|0,0,1,0,2,0>
0.069394
\n
|0,0,2,0,1,0>
0.061578
\n
|1,0,0,1,1,0>
0.059484
\n
|0,1,0,1,0,1>
0.053757
\n
|1,1,1,0,0,0>
0.042388
\n
|0,0,0,2,1,0>
0.039433
\n
|1,0,0,2,0,0>
0.039077
\n
|0,2,0,1,0,0>
0.037225
\n
|0,1,0,2,0,0>
0.031791
\n
|0,1,1,0,1,0>
0.028824
\n
|1,1,0,0,1,0>
0.028661
\n
|0,0,3,0,0,0>
0.028538
\n
|0,1,1,1,0,0>
0.024825
\n
|2,0,0,0,1,0>
0.024031
\n
|0,1,1,0,0,1>
0.022896
\n
|0,1,2,0,0,0>
0.021384
\n
|1,1,0,0,0,1>
0.021065
\n
|0,1,0,0,1,1>
0.020496
\n
|0,2,1,0,0,0>
0.020191
\n
|0,0,0,3,0,0>
0.020153
\n\n
"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"simulator_backend = pcvl.BackendFactory().get_backend(\"SLOS\")\n",
"simSLOS = simulator_backend(unitary)\n",
"\n",
"probsSLOS = pcvl.BSDistribution({state : prob for state, prob in simSLOS.allstateprob_iterator(input_state)})\n",
"pcvl.pdisplay(probs, max_v=20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also have to define the TVD function."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def tvd(probs1, probs2):\n",
" tvd = 0\n",
" for state, prob in probs1.items():\n",
" tvd += abs(prob - probs2[state])\n",
" return tvd\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we compute the TVD between the two simulations for different bond dimensions, going from 1 to 10."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": "",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEWCAYAAADRrhi8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvpklEQVR4nO3deXhV5bn38e+dhCRAwjwTZhBBRJCoKLai2ONQK9YZFXGoqBW1Pa0ebT2ttfb0rR2sCs5SnJFqWznW1tMqTqBAUAQRGQRkkBmZh5Bwv3+sFdzEnQmSvfbO/n2uK1f2mu/1rOFezxrN3REREUklGVEHICIiUlNKXiIiknKUvEREJOUoeYmISMpR8hIRkZSj5CUiIiknsuRlZm5mPQ9y2G+Y2YLajqka0+1tZrPNbJuZ3VTNYQ56PqNmZleY2btRxyE1Y2YTzOzuqOOoK4e471hmZqdW0O2A/Upsv2b2EzN7/OAirjKmruE8ZYXN/zCzUbU07grnqZbGP8/MhtbW+GqiyuQVzuwuM9se8zc2EcHFxHDAyuru77h770TGELoVmOLu+e5+f/mOZvammX2vticaroBlZb8jLI/Y5THezJ6KM9xRZrbHzFqY2Z1mtjdMvNvMbKGZjTWz9rUdbzjtpE98ZtYsLLs1MWVyW0z3CneSZnaWmc0Il8dGM3vWzApiumeb2e/NbGW4jJaZ2R8TMFuVCpeLm9m95doPD9tPCJvLdqjbY+K/rVz/s81sq5ltMLM3zKxbgmenVlW2X3H3/3H378HXk00dxHGGuz9ZVX/VSeK1ua+Md1Dk7ke4+5u1Mf6aqm7N6zvunhfzN6ZOo0peXYB5iZ5ouALmuXsecETYullMu8eAc82scblBRwKvuPumsPkFd88HWgDfBdoBs+oqgaWAe4E8oA/QFDgbWFzVQGZ2PvAc8EegFcEy2QO8a2bNw95uBwqBY4F8YCjwQa1Gf/A+Ay4st/MdBSyM02+zcB0bAfzMzE4Pd5hPAT8iKLduwDig9FADq6uEkI7qfVm6e6V/wDLg1Djtc4DNQL+Ydq2BXUCbsPkagp3BJmAy0CGmXwd6hr/fBL4X0+0K4N3w99thvzuA7cBFBDuClTH99wnHsZkguZwd020CwYb1d2AbMB3oUcn8nh2OY3M4zj5h+zcINs7dYRyHlRvuV+W6j42Zz+uAReE4xwEWM9xVwHzgS+A1oEsVy6NrOM6scu0XAJfHNGcCXwDDw+Y7gWfKDZMJfAT8roJpXQFMBcYCW4BPgWEx3ZsCTwCrgVXA3eE4+4TlUBqWxWaCHdxmICMc9jFgXcy4ngZ+UNl4q1NmVZV3ufn7GDinkrLev47GtDPgc+DWcu0zwvHdFTa/UjY/1fkD7gNWAFuBWcA3YrrdCUwiSBjbCNbPwpjuAwkS4zbgBWAicHcly/Rd4J/At8N2LYA1wG+BCRWtZ8BM4MfA+cDsGszbBOBh4F9hjG/FWWY3hMtsaTX3HTcBS4ANYdxl61UPgm11Y9jtWYIEHLs/ux34JFx//gTkht2GcuB+ZRnhvo+Y7QdYHsawPfw7KYzzyJhh2wA7gdZxyiMT+F0Y35Jw3veXNTH7Q6BnWF5bwv5fqGq/CPxXuDyfrmCeKpr/Kwj3u+W3AWA0sBcoDqf3v3HKKIfggO6L8O+PQE5s2RIc8Kwj2LavjJnOmWFM2wi2+R9XuV5VY8XbH1ycbuOBX8U03wD8M/x9SljYR4cz9QDwdrwdA5Ukr3g7kdgFAjQgWMl/AmSH090G9I7ZcDYSHAFnEazMEyuYn8PCleFb4XhvDcedHS/OOMN/rXsY+ytAM6AzsB44Pew2PBx/nzC2O4BpVSyPrsRPXj8F/h3TfFo4rQblN75yw90FTK9kR1cC/DAsj4sINqIWYfe/Ao8AjQk21hnAtZVsCMuBQeHvBQQbbp+YbgOrMd5Ky6yy8o4zf48TJIIrgV5xusdLXoeH7bvF6f8XwHvh7zvCefo+cCQVJNCYYS8DWobz9COCnU/ZTuVOgoOBMwl2fL8G3g+7ZRMk07JldD7BTqaq5HUJX+0Ivx+W993ESV4ECXsIwc54GNA9jOde4GQgr4p5m0CwTX6TYF9wH1/fvv9FkEQbUr19x5Sw/84ENcbYnf23wuFaE+zk/1huf/Yx0CkcfmpZWVH95LW/bGL6fRD4TUzzzYQ7+DjlcR3BgWBZDFOoOHk9T7BtZwC5wIlV7BdLgN+E89+wgnmqaP6voILkFbMc7y7XPbaM7gLeJ9hmWwPTgF+Wi+0ugvX0TIL1qXnYfTXhARvQHDi6snXKvfrJq+zouezvmrDbqcBnMf1OJTz6JzhyviemWx7BRtU1TqHsX1jxCrGChVSWvL5BsKFnxHR/HrgzpsAfj+l2JvBpBfP638CkmOYMgqOAofHijDP817qHsceucJOA28Lf/wCuLje9nVRS+6Li5NU5LN+CsPlZ4L6Y7ncSP3ldByyqYFpXEBxBxdYUZxCcjmxLcKqsYUy3EQTXBL+2DMN2TwP/SXC6cgFwTzj9/bWyaoy30jKrrLzjzF9DgoOeWWHZLQbOqGi9C9udGLbPrawsCZLMDQTbxJ6wHEdVtb3FjOtL4KiYZRd7YNIX2BX+/macZTSNqpNXQ2AtQS33fYLkFC95bQ5jmQ/cFDOewWHZridIZBOoIImF3SbGNOcR1Mo7xZTzKTHdq7PvOD2m+/eB1yuY9jnAhzHNy4DrYprPJNyHcWjJ6ziCgxULm4uACyuI6Y1yMfwHFSevp4BHCbfrcuOJt18sJmbdrGCeKpr/Kzi05PUZcGZMt9OAZTFx7CpXZuuAweHv5cC1QJPqbiPVveZ1jrs3i/l7LGw/BWhkZseZWVdgAMFRM0AHgiNCANx9O0ENqGM1p1ldHYAV7r4vpt3n5aazJub3ToKNoaJxxca8j+BUzqHGXNH0uwD3mdlmM9tMcOrBDmZ67r6c4CjzMjPLI9hov3YTRxwdw+lWZJWHa1foc4Jy6kJwBLU6Jv5HCI66KvIWwUr8zTDWNwlOuZwEvBOWd1XjrU6ZVWt5u/suDy7EDyKo9UwC/mxmLSqZhw3h/3jXCduXdXf3Uncf5+5DCGqBvwLGm1mfeCM1sx+b2Xwz2xLOV1OC62kVzVNueE2jA/GXUaXcfRfBqfQ7gJbuPrWCXlu5e3N37+MxNym5+/vufqG7tyY4gPwmQQ2hIitiht1OsNw6xOtO9fYdsf2XrZOYWVszm2hmq8xsK/AMB5ZjhcMeCnefTrBchprZ4QQ1wMkV9N4hTgwVuZVg/Z4R3tl3VRWhrHf33VX0U+vzHzpgucUZ90Z3L4lpjt02zyNIpJ+b2VtmdnxVEzukW+XdvZRggx8R/r3i7tvCzl8Q7GgACG8maElQkylvB9AoprldDcL4AuhkZrHz0rmC6VRnXLExG0H1urrj8qp7OcAKgtNhsQcGDd19Wg3HU+ZJglrReQTXDmZV1nNYZt8B3qmkt45hOZTpTFBOKwhqFK1iYm/i7mU3lMQri7cIdnRDw9/vEhzxnxQ2U43x1naZBcG6bwX+h+BUZWV3zS0gOHd/QWzLsCzPA16PM+5d7j6OoAbTt3x3M/sGwU7qQoLTKM0ITs9a+X7jWE38ZVQdZTddPFPN/uNy95nAX4B+lfTWqexHeHDVgmA92j+amN/V2Xd0ivldtk5CsAyd4PpTE4LTseXLsaJhq6ui7fzJcHojgRcrSSKr48QQf0Lua9z9GnfvQFAzebCKOwyrsw+qaP4P2A+bWfn9cFXjPmC5UYOydfeZ7j6c4CD1bwR5pVK18ZzXcwTXQi4Nf5d5HrjSzAaYWQ7BSjXd3ZfFGcdsgrvlGoUL5upy3dcSnGePp+yI51Yza2DBMwffIbhoXVOTgG+b2TAza0CwYe8hOA1THZXFGc/DwO1mdgSAmTU1swuqGKYyLxGsML8g2JDiMrOssAbwPMGBwh8qGWcb4KawbC8guNb0qruvBv4P+L2ZNTGzDDPrYWYnhcOtBQrMLLtsRO6+iODUwWXAW2HCWEuw038r7Keq8dZamZnZf5vZMeFt7bkE1yk2EySoMtlmllv2R7DN/Bi4w8wuCdu3I7h+1oTgOhBm9gMzG2pmDcPyHkVw1+GHcULJJ7gesB7IMrOfheOqjvfCYcuW0bkE13er4y2C60MPVLN/AMzsRDO7xszahM2HE9zo9H4lg50ZDpcN/JLgmt2KCvqtzr7jFjNrbmadCJbbC2H7fILLHFvMrCNwS5zx32BmBWEN+6cxw1bXemAfX9/WnyG4i/cyKj/rMYlgeRWEd6feVlGPZnaBffUIxpcECaTsLFNN9zdlKpr/j4AjwnLPJThVGquq6T1PsF20NrNWwM+oxoFRuP1damZN3X0vwU1L+6oarrrJ63/twOeKyk4NllWXdxBUD/8R0/7fBNeQXiI40ugBXFzB+O8lOFe7lmCn+2y57ncCT4anii6M7eDuxQTJ6gyCUzYPElx3+7Sa8xY7rgUEK94D4bi+Q/CYQHE1R3EfcL6ZfWlmX3sOLM70/kpwcXVieIrj43A+Doq77yAo7wK+XoYAF5nZdoKj+skEp2IGuXtlR0fTgV4E5fEr4Hx33xh2u5zghoGyO5de5KvTaW8Q3Ayxxsw2xIzvLYLTBytimo0DbyOvcLy1XGZOcLfVBoIjxG8R3IG3PaafeQQJt+zvSnd/geDo+ocEZfgJwTWkITFlsxP4PcHpvg0E17/Oc/clceJ4jeDuv4UEp1p2c+CpnYpnIFg3zyW4XrGJ4EDyL9Uc1t39df/qUYrq2kyQrOaG69M/CS4X3FPJMM8BPw9jHESwnVUUV3X2HS8TXKucTXD684mw/S8IbvTYEraPVxbPERwgLSG4TlOjB7rdfSfBtjA13CcNDtuvIFiPncrPZjxGsMw/CvuvbHkdA0wPy3kycHPMOnQnFewXqxB3/t19IcENFf8muPOz/HOaTwB9w+n9Lc547ya41jcHmBvOW3XLdiSwLNymryOoDFWq7OKiiEidsODB55XufkfUsdQ1MxsPfJEO8xq1+v0Qm4hIglhw09q5BM/dSR3Ti3lFRA6Rmf2S4BT2b919adTxpAOdNhQRkZSjmpeIiKSclL/m1apVK+/atWvUYYiIpJRZs2ZtCB8yT0kpn7y6du1KUVFR1GGIiKQUM6vyTSzJTKcNRUQk5Sh5iYhIylHyEhGRlKPkJSIiKUfJS0REUk7CkpeZjTezdWb2cQXdzczuN7PFZjbHzI5OVGwiIpJaElnzmgCcXkn3MwjeXt4LGA08lICYREQkBSXsOS93fzt8cWVFhgNPhV+Efd/MmplZ+/D7TrXuw+VfMmXBevJzsmick0VebtZXv8v+coP/2Vk6uyoikkyS6SHljhz4DaOVYbuvJS8zG01QO6Nz5+p+NPZAc1Zu4f7XF1Wr3+zMDPJys2ick0leToMwyWWSl9uAvJzMMNk1oHFOJvm5XyXA2N95YWJskKlEKCJyqJIpeVWbuz8KPApQWFh4UG8WHnVCVy4b3IWdxSVs31PCjj0lbNt94O8de4Lm7XtK2b5nL9t3f/V7w/Zilm3cGXTfXcKuvaXVmm5ug4wDklleThZ9OzThltN60yg7JReHiEjCJdPechXQKaa5IGxXZzIzjPzcBuTnNjjkcZWU7mNHcWn85Le7LAmG3cL/23eXsHX3XiZMW8bUxRt46LJB9GidVwtzJiJSvyVT8poMjDGzicBxwJa6ut5VF7IyM2jaMIOmDWueCN9dtIGbJn7I2Q+8y2/O789Z/TvUQYQiIvVHIm+Vfx54D+htZivN7Gozu87Mrgt7eRVYAiwGHgO+n6jYonZir1b8/aYT6d0unzHPfcidk+dRXLIv6rBERJJWyn+MsrCw0OvLW+WLS/bx63/M509Tl3F052aMu/Ro2jdtGHVYIlIPmdksdy+MOo6DpVvfkkh2VgY//84RjL1kIAvWbOPb97/Lu4s2RB2WiEjSUfJKQmf178DLY06kVV42I8dP5/7XF7FvX2rXkEVEapOSV5Lq2SaPv90whOFHdeAP/1rIVU/O5MsdxVGHJSKSFJS8klij7CzuvWgAd5/Tj2mLN3LWA+/y0YrNUYclIhI5Ja8kZ2ZcNrgLf77ueAAuePg9nn7/c1L9RhsRkUOh5JUijurUjFduPJETerbkv//2MT98YTY7i0uiDktEJBJKXimkeeNsxo86hh996zBe/ugLho+dyuJ126MOS0Qk4ZS8UkxGhnHjsF48fdVxbNxRzPCx7/LKnC+iDktEJKGUvFKU3sohIulMySuFtW/akImjj+fKIV2ZMG0ZFz/6Hqu37Io6LBGROqfkleLK3sox7pKj97+V451F66MOS0SkTil51RPf7t+eyTcGb+W4fPwMvZVDROo1Ja96pEdrvZVDRNKDklc9o7dyiEg6UPKqh/RWDhGp75S86jG9lUNE6islr3pOb+UQkfpIySsNxL6VY5PeyiEi9YCSVxo5sVcrXtFbOUSkHlDySjN6K4eI1AdKXmlIb+UQkVSn5JXGyr+VY9yUxbqdXkRSgpJXmit7K8fZR3Xgt68t4AcvzGb33tKowxIRqVRW1AFI9BplZ/HHiwZwWNt8fvvaAj7fuJNHLx9Em/zcqEMTEYlLNS8Bgrdy3HByTx6+bBAL1mzjnLFTmffFlqjDEhGJS8lLDnB6v3b8+brjceD8h97jtXlrog5JRORrlLzka/p1bMrLNwyhd7t8rn16lm7kEJGko+QlcbVpksvE0YP338jxn5M+0o0cIpI0dMOGVCi3QSb3XTyAXm3y+P2/FvL5xh08MrKQ1vk5UYcmImlONS+plFnwXsQHLz2aT1Zv5ZxxU5m/emvUYYlImlPykmo588j2/PnaEyjd55z30DT+9cnaqEMSkTSW0ORlZqeb2QIzW2xmt8Xp3tnMppjZh2Y2x8zOTGR8UrkjC5ry8pgh9GyTx+ini3j4rc90I4eIRCJhycvMMoFxwBlAX2CEmfUt19sdwCR3HwhcDDyYqPiketo2yeWF0cdz5pHt+X//+JRbXpzDnhLdyCEiiZXImtexwGJ3X+LuxcBEYHi5fhxoEv5uCuijU0moYXYmY0cM5Aen9uLFWSu57PHpbNy+J+qwRCSNJDJ5dQRWxDSvDNvFuhO4zMxWAq8CN8YbkZmNNrMiMytav15vQ4+CmfGDUw9j7CUDmbNyC8PHTWXBmm1RhyUiaSLZbtgYAUxw9wLgTOBpM/tajO7+qLsXunth69atEx6kfOWs/h2YdO3xFJfs49wHp/LGp7qRQ0TqXiKT1yqgU0xzQdgu1tXAJAB3fw/IBVolJDo5aEd1asbkMSfSrXVjrn6yiMffWaIbOUSkTiUyec0EeplZNzPLJrghY3K5fpYDwwDMrA9B8tJ5wRTQrmkuk649njP6tePuv8/nv16aQ3HJvqjDEpF6KmHJy91LgDHAa8B8grsK55nZXWZ2dtjbj4BrzOwj4HngCtchfMpolJ3F2BFHc9MpPZlUtJLLnpjOph3FUYclIvWQpXpuKCws9KKioqjDkHJenr2KW16cQ7smuTwxqpBebfOjDklEYpjZLHcvjDqOg5VsN2xIPTF8QEdeGD2YncWlnPvgNKYsWBd1SCJSjyh5SZ0Z2Lk5k8cMoVOLRlw9YSbj312qGzlEpFYoeUmd6tCsIS9efzzf6tuWu175hJ/89WP2lupGDhE5NEpeUucaZWfx0KWDuOHkHjw/Yzkjn5jOl7qRQ0QOgZKXJERGhnHLaYdz70VH8cHnmznnwaksXrc96rBEJEUpeUlCfXdgAc+PHsyOPSV898GpvL1Qj/GJSM0peUnCDerSnL/dMISOzRpy5YSZPDltWdQhiUiKUfKSSBQ0b8RL15/Ayb3b8PPJ87jjb3N1I4eIVJuSl0SmcU4Wj4wcxLUndeeZ95dzxZ9msGXn3qjDEpEUoOQlkcrMMG4/ow+/u+AoZizdxDkPTmXyR1+waO021cREpEJZUQcgAnD+oAK6tGzE9c/M4qbnPwSgQabRo3Uevdvlc1jbfHq3zad3u3w6NmtIRoZFHLGIREnvNpSksqeklM/W7WDh2m18umYbC9duY8GabazavGt/P42zM+nVNp/Dy5Ja+L91fk6EkYukllR/t6FqXpJUcrIy6duhCX07NDmg/bbde1m4dvv+ZLZgzTb+75O1TJz51ce5WzbO3p/MyhLaYW3zyM9tkOjZEJE6puQlKSE/twGDujRnUJfmB7TfsH3P/mS2cO02Fqzdxp+LVrCjuHR/Px2bNfzq1GO7PHq3bUKPNo3JycpM9GyISC1R8pKU1iovh1Y9cxjS86sPbu/b56zavOtrpx7fWbSevaXBafLMDKNry0Yc3q7J/qR2WNt8urRsTKaup4kkPSUvqXcyMoxOLRrRqUUjhvVpu7/93tJ9LNuwgwVrt7FwTZDY5n2xhVc/Xk3Zpd+crAx6tc3jiPZN+fFpvXUdTSRJKXlJ2miQmUGvtvnBhzH7f9V+V3Epi9dt59M1W8NTj9v5y4crAfjN+f0rGJuIREnJS9Jew+xMjixoypEFTfe3u3PyPJ55/3PGnNKTTi0aRRidiMSjh5RF4rj2pO5kmPHQW59FHYqIxKHkJRJH+6YNufCYAv5ctIIvYp4xE5HkoOQlUoHrTuqBOzys2pdI0lHyEqlAQfNGnD+ogIkzVrBmy+6owxGRGEpeIpX4/tCelLrzyNuqfYkkEyUvkUp0btmI7w7syHPTl7Num2pfIslCyUukCjec3JO9pft47O0lUYciIiElL5EqdGvVmOEDOvLM+8vZsH1P1OGICEpeItVyw8k92V1SyuPvLI06FBFByUukWnq2yeOs/h146r1lbNpRHHU4ImlPyUukmm48pSc7i0sZ/65qXyJRU/ISqabD2uZz5pHtmDBtGVt27o06HJG0ltDkZWanm9kCM1tsZrdV0M+FZvaJmc0zs+cSGZ9IVcac3Ivte0r40zTVvkSilLDkZWaZwDjgDKAvMMLM+pbrpxdwOzDE3Y8AfpCo+ESqo2+HJvxH37aMf3cpW3er9iUSlUTWvI4FFrv7EncvBiYCw8v1cw0wzt2/BHD3dQmMT6RabhrWi627S3hq2rKoQxFJW4lMXh2BFTHNK8N2sQ4DDjOzqWb2vpmdHm9EZjbazIrMrGj9+vV1FK5IfP06NmXY4W14/N2lbN9TEnU4Imkp2W7YyAJ6AUOBEcBjZtasfE/u/qi7F7p7YevWrRMboQhw47BebN65l6ff+zzqUETSUiKT1yqgU0xzQdgu1kpgsrvvdfelwEKCZCaSVAZ0asZJh7XmsXeWsLNYtS+RREtk8poJ9DKzbmaWDVwMTC7Xz98Ial2YWSuC04h6oZwkpZuG9WLTjmKefX951KGIpJ2EJS93LwHGAK8B84FJ7j7PzO4ys7PD3l4DNprZJ8AU4BZ335ioGEVqYlCX5pzYsxWPvL2EXcWlUYcjklbM3aOO4ZAUFhZ6UVFR1GFImpqxdBMXPvIePzurL1ed2C3qcESqzcxmuXth1HEcrGS7YUMkpRzbrQWDu7fg4bc+Y/de1b5EEkXJS+QQ3TSsF+u27WFS0YqqexaRWqHkJXKIju/ekmO6NuehNz9jT4lqXyKJoOQlcojMjJuG9WL1lt28OGtl1OGIpAUlL5FacGLPVgzs3IwHp3xGccm+qMMRqfeUvERqQVnta9XmXfz1Q9W+ROpajZKXmY0ysw/MbEf4V2Rml9dVcCKpZOhhrelf0JSxUxazt1S1L5G6VO3kZWajCD5R8iOgA8FLdW8FbjazkXUSnUgKMTNuPKUXKzbt4uXZX0Qdjki9VpOa1/XAd919irtvcffN7v4GcB5wQ92EJ5JaTu3Thj7tmzBuymJKVPsSqTM1SV5N3H1Z+ZZhuya1FZBIKjMzbh7Wk6UbdvDKnNVRhyNSb9Ukee06yG4iaeU/+rajd9t8HnhjEaX7Uvv1ayLJqibJq4+ZzYnzNxc4vK4CFEk1GRnGjcN68tn6Hbw6V7UvkbqQVYN++9RZFCL1zBn92tOzzSLGvrGYbx/ZnowMizokkXqlJjWvW4ECd/883l9dBSiSijIzjBtP6cmCtdv4v0/WRB2OSL1Tk+S1EPidmS0zs3vMbGBdBSVSH5zVvwPdWzXmvtcXk+qfHhJJNtVOXu5+n7sfD5wEbATGm9mnZvZzMzusziIUSVGZGcYNJ/dk/uqt/Hv+uqjDEalXavx6qPA04W/cfSAwAjiH4MvIIlLO8AEd6NyiEfe/vki1L5FaVOPkZWZZZvYdM3sW+AewADi31iMTqQeyMjMYc3JP5q7awpsL1kcdjki9UZPXQ33LzMYDK4FrgL8DPdz9Ynd/ua4CFEl13z26Ix2bNeQ+1b5Eak1Nal7/BKYBfdz9bHd/zt131FFcIvVGg8wMbji5J7NXbOadRRuiDkekXqhJ8prj7o+7+5d1Fo1IPXXeoI50aJqr2pdILanJQ8qtzew/K+ro7n+ohXhE6qWcrEyuH9qD/355Hu99tpETeraKOiSRlFaTmlcmkF/Jn4hU4oLCTrRtksN9ry+KOhSRlFeTmtdqd/9FnUUiUs/lNsjkupN68Iv//YT3l2xkcPeWUYckkrJqUvPSy9lEDtGIYzvTKi+HB95Q7UvkUNQkeQ2rsyhE0kRQ++rO1MUbKVq2KepwRFJWTV4PpS1NpBZcclxnWjbO5v43FkcdikjKqvEbNkTk0DTKzuKab3bn7YXr+XC5njwRORhKXiIRGDm4C80bNeAB1b5EDoqSl0gEGudk8b1vdOeNT9cxd+WWqMMRSTlKXiIRufz4LjTJzeJ+3XkoUmMJTV5mdrqZLTCzxWZ2WyX9nWdmbmaFiYxPJJHycxtw9Ynd+dcna/nki61RhyOSUhKWvMwsExgHnAH0BUaYWd84/eUDNwPTExWbSFSuGNKV/Jwsxk5R7UukJhJZ8zoWWOzuS9y9GJgIDI/T3y+B3wC7ExibSCSaNmzAFUO68urcNSxYsy3qcERSRiKTV0dgRUzzyrDdfmZ2NNDJ3f9e2YjMbLSZFZlZ0fr1+sCfpLarhnSjcXYmY6fozkOR6kqaGzbMLAP4A/Cjqvp190fdvdDdC1u3bl33wYnUoeaNs7n8hK68MucLFq/bHnU4IikhkclrFdApprkgbFcmH+gHvGlmy4DBwGTdtCHp4HsndiM3K5Nxqn2JVEsik9dMoJeZdTOzbOBiYHJZR3ff4u6t3L2ru3cF3gfOdveiBMYoEomWeTmMPL4LL89exdIN+kC5SFUSlrzcvQQYA7wGzAcmufs8M7vLzM5OVBwiyeqab3SnQWaGal8i1VCT73kdMnd/FXi1XLufVdDv0ETEJJIsWufncOlxXXjyvWXcdEovOrdsFHVIIkkraW7YEBG49qTuZGYYD76p2pdIZZS8RJJI2ya5jDimEy/OWsnKL3dGHY5I0lLyEkky1w3tQYYZD735WdShiCQtJS+RJNO+aUMuKCxgUtEKvti8K+pwRJKSkpdIErp+aA/c4ZG3VPsSiUfJSyQJFTRvxPmDCnh+5go+WrE56nBEko6Sl0iSuvnUXrTJz+HiR99nyqfrog5HJKkoeYkkqfZNG/KX759A99aN+d5TRUyauaLqgUTShJKXSBJrk5/LC9cezwk9WnLrS3O479+LcPeowxKJnJKXSJLLy8niiVHHcO7RHbn33wv5yV/nUlK6L+qwRCKV0NdDicjByc7K4PcXHEX7prmMm/IZ67bu4YFLBtIoW5uwpCfVvERShJlxy2mH88tz+jFlwTpGPDadjdv3RB2WSCSUvERSzMjBXXjoskF8unor5z00jc836hMqkn6UvERS0GlHtOO5a45j8669nPvgND0LJmlHyUskRQ3q0oKXrj+BhtmZehZM0o6Sl0gK69E6T8+CSVpS8hJJcXoWTNKRkpdIPaBnwSTd6CERkXpCz4JJOlHNS6QeKf8s2CV6FkzqKSUvkXqo7Fmw+XoWTOopJS+Reir2WbDzHprGnJWbow5JpNYoeYnUY2XPguU2CJ8FW6BnwaR+UPISqefKngXr1qox33uyiElFehZMUp+Sl0gaOOBZsBfncP/rehZMUpuSl0iaiH0W7A//0rNgktr0AIhIGtGzYFJfqOYlkmb0LJjUB0peImlKz4JJKlPyEkljehZMUlVCk5eZnW5mC8xssZndFqf7f5rZJ2Y2x8xeN7MuiYxPJB3pWTBJRQlLXmaWCYwDzgD6AiPMrG+53j4ECt29P/AicE+i4hNJZ3oWTFJNImtexwKL3X2JuxcDE4HhsT24+xR33xk2vg8UJDA+kbSmZ8EklSQyeXUEYg/nVobtKnI18I94HcxstJkVmVnR+vXrazFEkfSmZ8EkVSTlwx1mdhlQCJwUr7u7Pwo8ClBYWKhDQ5FapGfBJBUksua1CugU01wQtjuAmZ0K/BQ429318IlIBMo/C3bp49PZvqck6rBE9ktk8poJ9DKzbmaWDVwMTI7twcwGAo8QJC7d8iQSsZGDu/DgpUczZ+UWrn26iD0lpVGHJAIkMHm5ewkwBngNmA9Mcvd5ZnaXmZ0d9vZbIA/4s5nNNrPJFYxORBLk9H7tuee8/kxdvJEfvjCb0n06Uy/RS+hJbHd/FXi1XLufxfw+NZHxiEj1nDeogE07ivnVq/Np3uhj7j6nH2YWdViSxnQFVkSq5ZpvdmfjjmIefuszWuXl8MNvHRZ1SJLGlLxEpNr+6/TebNqxh/teX0TLvGwuP75r1CFJmlLyEpFqMzP+57tHsmnHXn4+eR7NGmVz9lEdog5L0pBezCsiNZKVmcHYSwZyTNcW/GjSbN5eqBcFSOIpeYlIjeU2yOTxUYX0bJPPdc/MYvaKzVGHJGlGyUtEDkqT3AY8edUxtMrL4co/zWDxuu1RhyRpRMlLRA5am/xcnr76WDIzMrj8iems3rIr6pAkTSh5icgh6dKyMROuPIZtu0sY+cQMvtxRHHVIkgaUvETkkPXr2JTHRhWyfNNOrnpyJjuL9R5EqVtKXiJSKwZ3b8kDIwby0YrNXPfMBxSX6FMqUneUvESk1px2RDt+fe6RvL1wPbe8+BH79B5EqSN6SFlEatVFx3Rm445i7vnnApo3yubn3+mr9yBKrVPyEpFad/1JPdi4vZgn3l1K6/wcbji5Z9QhST2j5CUitc7M+OmZfdi0o5jfvhbUwC45rnPUYUk9ouQlInUiI8O45/z+bN5ZzB1/m0vzRg0448j2UYcl9YRu2BCROtMgM4MHLx3EwM7NuXnibKZ9tiHqkKSeUPISkTrVMDuTJ0YV0rVVI0Y/NYuPV22JOiSpB5S8RKTONWuUzVNXHUfThg0YNX4GSzfsiDokSXFKXiKSEO2aBu9BdGDkE9NZu3V31CFJClPyEpGE6d46jwlXHsOXO4q5/IkZbNm5N+qQJEUpeYlIQvUvaMajlxeydMMOrn5yJruKS6MOSVKQkpeIJNyQnq2496IBzFr+JWOe+4C9pXoPotSMkpeIROLb/dtz1/B+vP7pOm57aS7ueg+iVJ8eUhaRyIwc3IVN24u5998LaZWXze1n9ok6JEkRSl4iEqmbhvVk0449PPL2Elo0zubak3pEHZKkACUvEYmUmfHz7xzBpp17+fU/PqV542wuLOwUdViS5JS8RCRyGRnG7y84is07i7n9L3Np3iibb/VtG3VYksR0w4aIJIXsrAwevmwQ/To2ZcxzHzBj6aaoQ5IkpuQlIkmjcU4Wf7riGAqaN+TqJ2cyf/XWqEOSJKXkJSJJpUXjbJ66+jjycrK4fPwMlm/cGXVIkoSUvEQk6XRs1pCnrjqWvaX7GDl+Ouu37Yk6JEkyCU1eZna6mS0ws8Vmdluc7jlm9kLYfbqZdU1kfCKSPHq1zWf8FcewbuseRo2fwdbdeg+ifCVhycvMMoFxwBlAX2CEmfUt19vVwJfu3hO4F/hNouITkeRzdOfmPHTZ0Sxcu43RTxWxe6/egyiBRN4qfyyw2N2XAJjZRGA48ElMP8OBO8PfLwJjzcxc740RSVtDe7fh9xcexc0TZ3Py794kL0dP+JS59qQenD+oIOowIpHItaAjsCKmeSVwXEX9uHuJmW0BWgIHfDvczEYDowE6d+5cV/GKSJIYPqAjAK/NWxNxJMmlWcMGUYcQmZQ8hHH3R4FHAQoLC1UrE0kDwwd03J/ERBJ5w8YqIPadLwVhu7j9mFkW0BTYmJDoREQkZSQyec0EeplZNzPLBi4GJpfrZzIwKvx9PvCGrneJiEh5CTttGF7DGgO8BmQC4919npndBRS5+2TgCeBpM1sMbCJIcCIiIgdI6DUvd38VeLVcu5/F/N4NXJDImEREJPXoDRsiIpJylLxERCTlKHmJiEjKUfISEZGUY6l+J7qZrQc+P8jBW1Hu7R1pTuVxIJXHV1QWB6oP5dHF3VtHHcTBSvnkdSjMrMjdC6OOI1moPA6k8viKyuJAKo/o6bShiIikHCUvERFJOemevB6NOoAko/I4kMrjKyqLA6k8IpbW17xERCQ1pXvNS0REUpCSl4iIpJy0TV5mdrqZLTCzxWZ2W9TxRMnMOpnZFDP7xMzmmdnNUccUNTPLNLMPzeyVqGOJmpk1M7MXzexTM5tvZsdHHVNUzOyH4TbysZk9b2a5UceUrtIyeZlZJjAOOAPoC4wws77RRhWpEuBH7t4XGAzckOblAXAzMD/qIJLEfcA/3f1w4CjStFzMrCNwE1Do7v0IPu2kzzZFJC2TF3AssNjdl7h7MTARGB5xTJFx99Xu/kH4exvBziltv7duZgXAt4HHo44lambWFPgmwbf2cPdid98caVDRygIahl96bwR8EXE8aStdk1dHYEVM80rSeGcdy8y6AgOB6RGHEqU/ArcC+yKOIxl0A9YDfwpPoz5uZo2jDioK7r4K+B2wHFgNbHH3/4s2qvSVrslL4jCzPOAl4AfuvjXqeKJgZmcB69x9VtSxJIks4GjgIXcfCOwA0vIasZk1JzhD0w3oADQ2s8uijSp9pWvyWgV0imkuCNulLTNrQJC4nnX3v0QdT4SGAGeb2TKC08mnmNkz0YYUqZXASncvq4m/SJDM0tGpwFJ3X+/ue4G/ACdEHFPaStfkNRPoZWbdzCyb4KLr5IhjioyZGcE1jfnu/oeo44mSu9/u7gXu3pVgvXjD3dP26Nrd1wArzKx32GoY8EmEIUVpOTDYzBqF28ww0vTmlWSQFXUAUXD3EjMbA7xGcMfQeHefF3FYURoCjATmmtnssN1P3P3V6EKSJHIj8Gx4oLcEuDLieCLh7tPN7EXgA4I7dD9Er4mKjF4PJSIiKSddTxuKiEgKU/ISEZGUo+QlIiIpR8lLRERSjpKXiIikHCUvSWlmVmpms83sIzP7wMxq5aFRMxtanTfKm9mbZlYY/n7VzJrVxvRrwszuMrNTEz1dkSil5XNeUq/scvcBAGZ2GvBr4KQoAnH3MyOa7s+imK5IlFTzkvqkCfAlBG8NMbPfht9dmmtmF4Xth4a1pbLvUz0bvi2h7Btvn5rZB8C58SZgZg3NbGL4Xau/Ag1jui0zs1Zm1jUczwQzWxhO41Qzm2pmi8zs2LD/xmY23sxmhC+9HR62v8LM/mJm/wz7vydsnxmOs2yefhi2n2Bm54e/h4XjmhuOOycmtl+EtdO5ZnZ4nSwBkQRRzUtSXcPwrSC5QHvglLD9ucAAgu9PtQJmmtnbYbeBwBEEn7OYCgwxsyLgsXD4xcALFUzvemCnu/cxs/4Eb1uIpydwAXAVwevILgFOBM4GfgKcA/yU4PVTV4WnG2eY2b/D4QeEce4BFpjZA0AboGP4LSnKn6IMP4w4ARjm7gvN7Kkw3j+GvWxw96PN7PvAj4HvVRC7SNJTzUtS3S53HxB+KPF04KmwJnUi8Ly7l7r7WuAt4JhwmBnuvtLd9wGzga7A4QQvXV3kwWtnKnoZ7zfLurn7HGBOBf0tdfe54TTmAa+H450bTg/gP4DbwuT7JkEC7hx2e93dt7j7boJ3CXYheDVTdzN7wMxOB8q/+b93ON2FYfOTYbxlyl64PCsmBpGUpOQl9Ya7v0dQy2pdRa97Yn6XUjdnIGKnsS+meV/M9Aw4L0y+A9y9s7vPjzN8KZDl7l8S1CTfBK6j5h/LLBtnXc2zSMIoeUm9EV7HyQQ2Au8AF4XXiVoT1EBmVDL4p0BXM+sRNo+ooL+3CU4BYmb9gP6HEPJrwI0x19wGVtazmbUCMtz9JeAOvv5pkgUE89AzbB5JUOMUqXd09CWpruyaFwQ1mVHuXhreTHE88BHgwK3uvqaiGxXcfbeZjQb+bmY7CZJffpxeHyL4qvB8gs9hHMpHK39JcD1qjpllAEuBsyrpv2M47bKDztvjzMOVwJ8t+Ez9TODhQ4hPJGnprfIiIpJydNpQRERSjpKXiIikHCUvERFJOUpeIiKScpS8REQk5Sh5iYhIylHyEhGRlPP/AeiFCeevQTAPAAAAAElFTkSuQmCC\n"
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"TVD = []\n",
"for chi in range(1,11):\n",
" simMPS.set_cutoff(chi)\n",
" probsMPS = pcvl.BSDistribution({state : prob for state, prob in simMPS.allstateprob_iterator(input_state)})\n",
" TVD += [tvd(probsMPS, probsSLOS)]\n",
"\n",
"plt.plot(TVD)\n",
"plt.xlabel('Bond dimension')\n",
"plt.ylabel('TVD')\n",
"plt.title('Evolution of the TVD between SLOS and MPS probability distributions');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the TVD decreases as the size of the matrices increases, untill reaching 0 for a bond dimension of 7."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"> [1] Ulrich Schollwöck. The density-matrix renormalization group in the age of matrix product states. [Annals of Physics](https://doi.org/10.1016/j.aop.2010.09.012), 326(1):96–192, jan 2011.\n",
"\n",
"> [2] Changhun Oh, Kyungjoo Noh, Bill Fefferman, and Liang Jiang. Classical simulation of lossy boson sampling using matrix product operators. [Physical Review A](https://doi.org/10.1103/PhysRevA.104.022407), 104(2), aug 2021.\n",
"\n",
"> [3] Hui Wang, et al. Boson Sampling with 20 Input Photons and a 60-Mode Interferometer in a $10^{14}$-Dimensional Hilbert Space. [Physical Review Letters](https://link.aps.org/doi/10.1103/PhysRevLett.123.250503), 123(25):250503, December 2019. Publisher: American Physical Society."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.4 64-bit",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.4"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "256f899c2a79c403df4e9a617fa483f1c8ef03431534114d54980d272fd739d5"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}